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We present an iterative sampling method which delivers upper and lower bounding processes 
for the Brownian path. We develop such processes with particular emphasis on being able to 
' unbiasedly simulate them on a personal computer. The dominating processes converge almost 

surely in the supremum and Li norms. In particular, the rate of converge in Li is of the order 
0{1C~^^^), K. denoting the computing cost. The a.s. enfolding of the Brownian path can be 
exploited in Monte Carlo applications involving Brownian paths whence our algorithm (termed 
the e-strong algorithm) can deliver unbiased Monte-Carlo estimators over path expectations, 
overcoming discretisation errors characterising standard approaches. We will show analytical 
results from applications of the e-strong algorithm for estimating expectations arising in option 
pricing. We will also illustrate that individual steps of the algorithm can be of separate interest, 
^ ' giving new simulation methods for interesting Brownian distributions. 
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1. Introduction 



Brownian motion (BM) is an object of paramount significance in stochastic modelling. 
Starting from its original mathematical formulation by [2], its properties are still under 
meticulous investigation by contemporary researchers. Relevant to the purposes of this 
^ ' paper, considerable work has focused on various constructions and representations of BM 

paths. Leaving aside the simple finite-dimensional Gaussian structure of BM, researchers 
have often been interested on more complex functionals. Hitting times, extremes, local 
times, reflections and other characteristics of BM have been investigated (for a general 
exposition see [18]). For simulation purposes, many of the relevant distributions are easy 
to sample from on a computer ([10]). Several conditioned constructions of BM arc also 
known relating BM with the Besscl process, the Raylcigh distribution and other stochastic 
objects (see e.g. [3]). 

This paper presents a contribution of our own at simulation methods for Brownian 
dynamics. We develop an iterative sampling algorithm, the e-strong algorithm, which sim- 
ulates upper and lower paths enveloping a.s. the Brownian path. To meet this objective, 
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we collect a number of characterisations and combine them in a way that they can de- 
liver simple sampling methods implementable on a personal computer. We will show that 
after C'(/C)-computational effort, the dominating process have Li-distance of 0{IC~^^^). 
This a.s. enfolding of the Brownian path can be exploited in Monte Carlo applications 
involving Brownian motion integrals, minima, maxima or hitting times; in such scenaria, 
the e-strong algorithm can deliver unbiased Monte-Carlo estimators over Brownian ex- 
pectations, overcoming discretization errors characterising standard approaches (for the 
latter approaches, see for instance the exposition in [12] in the context of applications in 
finance). 

We will show applications of the algorithm and experimentally compare the required 
computing resources against typical alternatives employed in the literature involving 
Euler approximation. Our examples will involve a collection of double-barrier option 
pricing problems in a Black and Scholes framework arising in finance. Also, we will 
demonstrate that individual steps of the algorithm can be of separate interest, giving 
new simulation methods for interesting Brownian distributions. 

The e-strong algorithm delivers a pair of dominating processes, denoted by X^{n) = 
u G [0, 1]} and X'^{n) = {X^(n); u G [0, 1]}, that can be simulated on a personal 
computer without any discretisation error, with the property: 

Xiin) < Xiin + 1) < X„ < Xlin + 1) < Xl{n) (1.1) 

for all instances u G [0, 1]; here, X is the Brownian path. The two dominating processes 
will converge in the limit: 

w.p.l , lim sup \xl{n) - Xi{n)\ ^ . (1.2) 
"^°°ue[o,i] 

The algorithm builds on the notion of the intersection layer, a collective information, 
containing the starting and ending points of a Brownian path together with information 
about its extrema. A number of operations (bisection, refinement, see main text) can 
be applied on this information, explicitly on a computer, allowing the sampler to iterate 
itself to get closer to X. 

We should note here that the methods described in this paper will be relevant also for 
non- linear Stochastic Differential Equations (SDEs). Recent developments in the simula- 
tion of SDEs under the framework of the so-called 'Exact Algorithm' (see [8, 7, 5, 13, 6, 9]) 
build upon the result that, conditionally on a collection of randomly sampled points, the 
path of the SDE is made of independent Brownian paths. Once this collection of points 
is sampled, the methodology of this paper can then be applied separately on each of the 
constituent Brownian sub-paths. 

The structure of the paper is as following. In Section 2 we present the notion of 
the intersection layer which will be critical for our methods. In Section 3 wc present the 
individual steps forming the e-strong algorithm; they will require original simulation tech- 
niques for some Brownian distributions. Once we identify in Section 4 the (-function, an 
alternating monotone scries at the core of Brownian dynamics, we exploit its structure in 
Section 5 to analytically develop these new sampling methods. In Section 6 we apply the 
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e-strong algorithm to unbiasedly estimate some path expectations arising when pricing 
options in finance. We will contrast the computational cost of the algorithm with Euler 
approximation alternatives to get a better understanding of its practical competitiveness. 
In Section 7 we sketch some other potential applications of the e-strong algorithm. We 
finish with some discussion and conclusions in Section 8. 

2. Intersection Layer and Operations 

We will, in general, write paths as X = {X^; u S [s, <]} for s <t. A Brownian bridge on 
[s, t\ is a Brownian motion conditioned to start at Xg and end at Xt, for some pre-specified 
Xs, Xt] its finite-dimensional dynamics are easily derivable following this interpretation 
(see for instance [18]). 

Instrumental in our considerations is the notion of (what we call) the intersection 
layer. Consider a Brownian bridge X on [s,i\. Let rris.t, Mg.t be the extrema of X: 

rus^t = ird{Xu\u G [s,i]} , Ms.t = sup{X„;u G [s,t]} . 

The e-strong algorithm will require some information on both m^^t and M^^t- We will 
identify intervals: 

such that: 

We will write simply to, Af, , U'^, L^, ignoring the s, t-subscripted versions when 
the time interval under consideration is clearly implied by the context. The intersection 
layer idea refers to the collective information 

^sj ~ {Xs, Xt, Cs,t,l^s,t} , (2-1) 

that is the starting and ending points of the bridge together with intervals that contain 
its maximum and minimum. Fig. 1(a) presents a graphical illustration of the intersection 
layer: the extrema of an underlying Brownian bridge lie in the shaded rectangles. We will 
look now at two simple operations on the information I,,* which nonetheless will be the 
building blocks of the complete e-strong algorithm described in the next section. 

2.1. Refining the Information Is,t 

During the iterations at the execution of the e-strong algorithm, for each piece of infor- 
mation Is^t we will need to control the width of the layers Cg^t-, l^s,t relatively to the size 
t — s oi the time interval to ensure convergence of the bounding paths enveloping the 
underlying Brownian path, Thus, the refinement of the information Ig^t corresponds to 
a procedure that updates Xs,t by halving the allowed width for the minimum to or the 
maximum M of the path, thereby correspondingly updating the layers Cg^t or Us.t- 
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Bisect(Ts,t ): 

1. Set t* = {t + s)/2. Simulate Xf given Xs,t- Set (7-1- = ^7■^ V Xf, , = L'^ A Xf . 
2a. Decide ifUs,f = [Xs V Xt* , or [U^, C/t], 
2b. Decide ifUf.t = [Xf V Xt,U^ or [U^,U^]. 
2c. Decide if Cs,f = [L^,L^ or [L^^Xs hXf]. 
2d. Decide if Cf ,t = [L^^L^ or [L^.Xf hXt\. 

3. Return III t * \J If t. 



Table 1. The procedure for bisecting the information Is.t- It returns the intersection layers (• and 
If .t with refined information about the underlying path (compared to Is,t). 

More analytically, refinement of Xs t corresponds to deciding whether the minimum m 
on [s,i], already known to be in lies in [L^, [L^ + L^)/2] or [{L^ + L'^)/2,L\ 

i.e. whether Cs,t is equal to [L^, {L^ + i^)/2] or [{L^ + L'^)/2, L^]; the apparent analogue 
of such a consideration applies for the maximum M. The analytical method of sampling 
the relevant binary random variables for carrying out this procedure will be described in 
Section 5. 

2.2. Bisecting the Information Is,t 

This is a more involved operation on Is,t, and involves bisecting Ts^t into the more 
analytical information Is,f Vlf^t for some intermediate time instance t* £ (t,s). In 
particular, wc will be selecting t* = {t + s)/2 within the e-strong algorithm. The method 
begins by sampling the middle point Xf conditionally on Is,t, and then appropriately 
sampling the layers for the two pieces of information 2s. t*, 2f,t- The practicalities of 
implementing the second part of the method will depend on whether Xf falls within a 
layer of Ts,t or not, thus we present the bisection operation in more detail in Table 1. 

Note that if Xf > U'^ the two upper layers (for 2s.f and If,t) will be directly set 
to [Xf,U^], and we will have to simulate extra randomness about the underlying path 
only to determine the lower layers. Correspondingly, if Xf < the two lower layers 
will immediately be set to [L^,Xf]. In the scenario when < Xf < we will have 
to simulate extra randomness to determine all four layers. We describe in Section 5 
the algorithms for sampling Xf and determining the layers. Fig.l shows a graphical 
illustration of the bisection procedure. 
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(a) 




(b) 



Figure 1. Top panel (a): the intersection layer information Xa,t for a Brownian path. The underlying 
trajectory starts at Xs and finishes at Xt with its extrcma found in the shaded areas. Bottom panel 
(b): the bisection of Is,t into Is^t* and Xf^t- The algorithm simulates Xf and then decides that the 
extrema for each of the intervals {s,t*] and {t* ,t\ are in the shaded areas, i.e. 11^ t* = l^f ^U'^], Uf ,t = 
[U^,U^], Cs,t' = [L^,L^] and Cf,t = lL^,Xf]. The algorithm outputs the upgraded information 
Is,t' = {^s,Xt*,Cs^t'Ms,t*] andXf^t = {Xf , Xt, Cf ,tMf ,t}- 
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E-strong(Xo, Xi, n): 

1. Initialize Uo^i, £0,1, set Io_i = {Xq, Xi,Uo,i, Co.i}- Set V = {Io,i} ind i = l. 

2. For each of the 2^~^ intersection layers in V , say Xs.t, do the following: 

i. Bisect the information Is,t into X^ -ff, Xf t, where t* = {t + s)/2. 
a. Refine X^ (», It*,f until the width of their layers is not greater than <^{t — s)/2. 

3. Collect the updated information, V = Vj=i ^{j-i)2-» j2-' ■ 

4- If i < n set i = i + 1 and return to Step 2; otherwise return V . 



Table 2. The e-strong algorithm. It itcratively unveils extra information about the underlying path. It 
outputs the collection of intersection layers V = Vj=i ^(j — 1)2-" j2-" • 

3. e-Strong Simulation of Brownian Path 

We introduce an iterative simulation algorithm with input a Brownian bridge X on 
the domain [0, 1] and output, after n iterations, upper and lower dominating processes 
X^{n) = {X^{n)] u g [0, 1]} and X^{n) = {Xl{n)]u e [0, 1]} satisfying the monotonicity 
and limiting requirements (1.1) and (1.2) respectively. Note that X here is a continuous 
time Brownian bridge path, thus an infinite-dimensional random variable. However, the 
bounding processes will be piece-wise constant, thus inherently finite-dimensional. One 
will be able to realise complete sample paths of X-'-(n) or X'^{'n) on a computer without 
retreating to any sort of discretization or approximation errors (apart from those due to 
finite computing accuracy). 

3.1. e-Strong Algorithm 

Given some initial intersection layer information 2o,ij the algorithm will naturally set 
X2(0) = Uq ^ and X^{0) = Lq i for all instances u G [0, 1]. It will then iteratively bisect 
the acquired intersection layers, as described in Section 2.2, to obtain more information 
about the underlying sample path on finer time intervals. To ensure convergence of the 
discrepancy X^(n) — X^{n) the algorithm will sometimes refine the information on some 
intersection layers, as described in Section 2.1, to reduce the uncertainty for the extrema. 
We give the pseudocode about the algorithm in Table 2. 

Utilising the information the e-strong algorithm returns, we define the dominating 
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processes as follows: 

2" 

^lin) = 5I^(l-l)2-",»2-" •I«G((»-l)2-".»2-"] , 

(3-1) 

^ii'f^) X]^('i-l)2-",i2-" ■ I«e((i-l)2-",i2-"] ■ 
1=1 

The square-root rate at Step 2.ii of the algorithm in Table 2 is to guarantee convergence 
of the dominating paths with minimal computing cost: it provides the correct distribu- 
tion of effort between time-interval and extrema-intcrval bisections. To understand this, 
note that the range of a Brownian motion (or a Brownian bridge) on [0, 2~"] scales as 
0(2~"/^); see for instance [18]. Thus, had we used the actual Brownian minima and 
maxima to define dominating processes for the Brownian path in the way of (3.1) the 
rate of convergence would have been 0(2~"/^); we cannot exceed such a rate, but we 
can preserve it if our extrema are not further than 0(2~"/^) from the actual ones. This 
intuitive statement will be made rigorous in the sequel, when an explicit result on the 
rate of convergence of the dominating processes in Li-norm is given. 

Fig. 2 shows successive steps of the e-strong algorithm as implemented on a computer. 
For each n, the red lines show the interval where the maxima and the minima are located: 
this information is available for all 2" sub-intervals bisecting the initial time interval [0, 1]. 
The black line corresponds to the linear interpolation of successively unveiled positions 
of the underlying Brownian path. The last graph (f) corresponds to n — 12; in this case, 
we have zoomed on a particular subinterval of [0, 1] to be able to visualise the difference 
between the bounding paths and the underlying Brownian one. 

3.2. Convergence Properties 

Almost sure convergence of the dominating paths follows directly from the continuity of 
the Brownian path X. The analytical proof is given in the following proposition. 

Proposition 3.1. Consider the continuous-time processes X'^{n), X^{n) defined in 
(3.1). Then, the convergence in supremum norm in (1.2) will hold in the limit n ^ oo. 

Proof. For a Brownian bridge X on [0, 1] we consider: 

£)„ := sup (M(i_l)2-'«,j2-" - '7l(,_i)2-",i2-") ■ 
l<i<2" 

Uniform continuity implies that, with probability 1: 

lim Dn = Q . 

Now, we have that: 

sup \Xl{n) - Xiin)\ < A. + 2 • 2-"/^ ^ , 
ue[o,i] 
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(a) n = 



(b) n = 1 




(c) n : 



(d) n : 




(e) n = 4 



(f) n = 12 



Figure 2. The e-strong algorithm as applied on a personal computer. For each step n, the red lines show 
the allowed interval for the minima and the maxima: this information is separately available for all 2" 
time sub-intervals partitioning [0, 1]. Note that the last graph corresponds to n = 12, with the subplot 
in its frame corresponding to a zooming on the position of the paths on the time interval [0.424, 0.434] . 
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where we have used the fact that Step 2.ii of the e-strong algorithm guarantees that 

^ {1-1)2-", i2~" - '^(»-l)2-",i2-" - ^ 



□ 



A more involved result can give the rate of convergence of the dominating processes and 
will be of practical significance for the efficiency of Monte-Carlo methods based on the 
e-strong algorithm. 

Proposition 3.2. Consider the Li-distance: 

\X\n) - = / \Xl{n) - Xi{n)\du . 

Jo 

Then: 

2"/2 X E [ \X^{n) - X^{n)\i] = . 
Proof. Wc proceed as follows: 

\xHn) - xHn)h = i:(C/(ii)2-..2-.. - 4-i)2-".2-0 ■ 

2" 

<E(^%-i)2-",.2-" -"^(«-i)2-N.2-+2-2-"/2) .2-" , (3.2) 

i=l 

the incciuality being a direct consequence of Step 2.ii of the e-strong algorithm in Table 2. 
Consider now the path from X(i_i)2-" to Xj2-n- Let Z he a Brownian bridge from 
= to ^2-11 = 0; we denote by and its maximum and minimum respectively. 
Conditionally on ^(i_i)2-" and X^2-", a known property of the Brownian bridge implies 
(see e.g. [14]) that: 

Xt+{i-i)2-" = Zt + (1 - -^(.i-i)^'" ^ 2^ -^t^'"' ie[0,2 "] , 

in the sense that the processes on the two sides of the above equation have the same 
distribution. It is now clear that: 

A'/(i-l)2-",i2-" - "l(i_i)2-",i2-" < - X(j_i)2-"| + (^-4 - m^) . 

So, taking expectations at (3.2), we get: 

2" 

E [ \X^n) - X\n)\,] <E[M,- ] + 2 • 2^"/^ -f ^E |X,2-n - X(,_i)2-,. | 2"" 
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The finite-dimensional distributions of the initial Brownian bridge from to Xi imply 
that: 

^.2- - ^(.-1)2-" ~ ^((^1 - ^o)2-", 2-«(l - 2-")) , 
which gives directly that: 

It remains to show that ¥j[Mz — niz] = 0(2^"/^) to complete the proof. Now, self- 
similarity of Brownian motion implies that: 

where Z is a Brownian bridge from Zq = to Zi = 0. Let M^, fhz be the maximum and 
minimum of Z. Due to the self-similarity we have 

M, - = 2-"/2(A4 - rh,) . 

Since Mz — rhz in a random variable of finite expectation (see e.g. [14]), we obtain directly 
that E[Mz -niz]^ 0(2""/^) which completes the proof. 

□ 



4. The (^-Function 

We have yet to present the sampling methods employed when refining or bisecting an 
intersection layer during the execution of the e-strong algorithm, thus constituting the 
building blocks of our algorithm. All probabilities involved in these methods can be 
expressed in terms of a hitting probability of the Brownian path. We denote by 

the probability law of a Brownian bridge from Xq = x to Xi = y. Let ^{L, U;l,x,y)^ with 
L < U, he the probability that the Brownian bridge escapes the interval [L, U]. That is: 

CiL,U;l,x,y)^W^'^''^y'> [mnj < L or Ah.i > U] . 

We also define: 

jiL,U;l,x,y)^l-aL,U;l,x,y) . (4.1) 

These probabilities can be calculated analytically in terms of an infinite series. The result 
is based on a partition of Brownian paths w.r.t. to a trace they leave on two bounding 
lines and can be attributed back to [11]; for more recent references see [1, 17, 9]. We 
define for j > I, 

a^{x,y,5,C) = exp |-y [5j + C - x] [Jj + C - 2;]| , 

f. (x, y, 5) = cxp I - ^ [ + 5{x - y) ] I . (4.2) 
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Then, Theorem 3 of [17] yields 

C(L,C/;^x,2/) = | L<x,y<U^ 

' ' ' \ 1 , otherwise , ^ ' 

where 

(Tj = CT. {x, y,U - L,L) + {-x, -y, U - L, -U) , 

Tj = (x, y,U-L)+f^ i~x, ~y, U ~ L) . (4.4) 

The infinite series in (4.3) exhibits a nionotonicity property which wih be exploited by 
our simulation algorithms. We consider the sequence {Sn}, with Sn — Sn{L,U;l,x,y), 
defined as: 

52,1-1 = ^^{<^3 - 'Tj) + O'ri, 52,1 = S'2n-1 - Tn . (4.5) 

i=i 

when L < x,y < U, otherwise 5„ = 1. Then: 

< S2n < S2n+2 < C < S2n+1 < S2n-1 (4.6) 

for all n > 1; for a proof see [9] or [6]. 
4.1. (^-Derived Events 

We can combine (^-probabilities to calculate other conditional probabilities arising in the 
context of the e-strong algorithm. We begin with the following definition: 

P{L^,L^,U^,U^]l,x,y) := w(''-^''«) [L^ < mo,; < L^, < Mq^i <U^] . 

Now, we have the set equality: 

{L^ < mo,/ < L\ < Mo,; <U^ = {L^ < nioj, M^j < U^} - (4.7) 

{L^ < moj, Moj < U^} U {L^ < moj, Mqj < U^} . 

Thus, taking probabilities and recalling the definition of 7 in (4.3), we find that: 

I3{L\ L\ U\ U^; /, X, y) - ^{L\ [/t; I, x, y) - j{L\ f/t; /, x, y) 

~^iL\U^;l,x,y)+j{L\U^;l,x,y) . (4.8) 

Before the next event, we enrich the notation for the Brownian bridge measure. We define 
(for 0<q<l): 

We set r ^ I — q. Consider now the conditional probability: 

p{L\L\u\U^;q,r,x,w,y)=wf^^j;^ [L^ < mo,/ < L\ < Mo,/ <U^]. 
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Using again the set equality (4.7), and taking probabilities under W|^'^^ , we obtain: 

p{L^, L^, U^, U^; q, r, x, w, y) = 7172 - 7374 - 7576 + 7778 (4.9) 
where we have defined: 

71 = l{L^, U^; q, X, w) , 72 = 7(L^, U^; r, w, y) , 73 = l{L'^, U^; Q, x, w) , 
74 = liL"^, U^; r, w, y) , 75 = l{L^. U^] 9, x, w) , -fe = l{L^, U^] r, w, y) , 
77 = liL'^, U^; q, X, w) , 78 = j{L'^, U^; r, w, y) . 

Note that the product terms arise due to the independency of the Brownian bridges on 
[0, q] and [q, I]. We will be using these expressions for /?(•; •) and p(-; •) in the sequel. 

4.2. Simulation of ^-Derived Events 

We will need to be able to decide whether events of probability ( have occurred or not. 
In a simulation context, this corresponds to determining the value of the binary variable 
Ifi<C ^ ^ Un [0, 1]. With (4.6) in mind, we define: 

J = inf{n > 1 : 71 odd, 5„ < i? or n even, S'„ > R} . 

Due to the alternating monotonicity property (4.6) of Sn- 

l-i?<{^ — I J is even ■ 

Thus, we need a.s. finite number of J computations to evaluate 1r<(. Note that Sn 
converges to its limit exponentially fast, so J will be of small expectation; one can easily 
verify that all its moments are finite. Such an approach was also followed in [6] . 

In a more general context, we will also be required to decide if events of probabil- 
ity /3(-;-) or /o(-; •) have taken place or not; we will in fact be considering even more 
complex events related with the ^-function. In the most encompassing scenario, when 
executing our sampling methods, we will be required to compare a given real number R 
with Z(Ci, C2, • • ■ , Cm) for some given function Z, with the different Q's corresponding to 
different choices of the arguments I, x, y, L, U for (■(■;■)■ Using the monotonicity property 
(4.6) we will be able to develop corresponding alternating sequences such that: 

^ ^2n+2 ^^iClX2, • • • : Cm) < < ! 

hm 5,f = Z(Ci,C2,...,Cm) , 

n— >oo 

and proceed as above. Analytically, we will determine the value of the comparison binary 
indicator I/f<z(Ci,c-3,- --Cm) follows: 

Calculate until the first n such that either n is odd and < R ( whence return ) 
or n is even and > R (whence return 1). 
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We will now describe analytically all simulation algorithms employed at the development 
of the e-strong algorithm presented in Table 2. In particular, one has to develop sampling 
methods to carry out the refinement and bisection (see Section 2) of the intersection layer 
Ts,t- To simplify the presentation, when conditioning on Xs, Xt* or Xt we will make the 
correspondence: 

X ^ Xs, w ^ Xt. , y ^ Xt , 
I ^t- s, q = t* - s, r = t-t* . 

5.1. Bisection of Xs,t' Sampling the Middle Point Xt*. 

Bisection oi Is,t = {Xs, Xt, Cs,t,l^s,t}, with Cs,t = [L^,L'^], Us^t = [U^,U^], begins by 
sampling a point of the Brownian bridge conditionally on the collected information about 
its minimum and maximum; this is Step 1 of Table 1. Such a conditional distribution is 
analytically tractable via Bayes' theorem. 

Proposition 5.1. The distribution W[Xt* | Is,*], with t* G [s,t\, has probability den- 
sity: 

f{w) oc p{L^, L^, U^, U^; q, r, x, w, y) x tt{w) 
where •) is defined in (4-9) and 

4«^)=cxp{-i(«;-(rx+fy))V(f)} . 

Proof. The function 7r(w) corresponds to the prior (unnormalised) density for the middle 
point Xt I Xs , Xt which is easily found to be normally distributed with mean and variance 
as implied by the expression for tt{w). So, following the definition of p{-\ ■) in (4.9), the 
stated result is an application of Bayes' theorem. □ 

We will develop a method for sampling from f(w). It is easy to construct an alternating 
series bounding f{w). Let: 

C^^l-l^ , 1 < « < 8, 

for the eight 7-functions appearing at the definition of p in (4.9). Let {>S'i,„}„>i be the 
alternating series (4.6) for d, for each 1 < i < 8; that is: 

< Si,2n < 'S'i,2n+2 < Ci ^ Si^2n+1 < <S'i,2n-l , (5-1) 

with lim„_j.oo Si^n = 0- Consider the sequence {S^} defined as follows: 

Sn = (1 ^ -^Ln+l — S2,n+1 + 'S'l,ri'5'2,n) ^ (1 ^ S^^i — 5*4, „ + Ss^n+lS^^n+l) 

— (1 — S^^n — S5^n + •S's.n+l'S'e.n+l) + (1 ~ S^^n+l — S's.n+l + 'S'y.n'S's.n) • (5-2) 
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Due to (5.1), one can easily verify that {S',f } is an alternating sequence for p(-; •), in the 
sense that: 

SL < < P{L\ L\ U\ U\ q, r, X, w, y) < S^^, < Sl_, (5.3) 

with lini„^oo = p(L^, L^, U^, If^; q, r, x, w, y). 

We exploit this structure to build a rejection sampler to draw from the density f{w) 
in Proposition 5.1. We will use proposals from: 

/2«+i(w) = S^^^^{w) X ■n{w) , 

where we have emphasized the dependence of S'^,j_|_i on the argument w. Note that the 
domain of both /(w), f2n+i{w) is [L^^W]. Now, we will illustrate that S2n+i{w) has 
a concrete structure that we will exploit for our sampler. Consider the first of the four 
terms forming up 5*^1+1 from (5.2): 

1 — <5'l,2n+2 — <5'2,2n+2 + S'l,2n+lS'2,2n+l • (5-4) 

Following the analytical definition of the alternating sequences in equations (4.2), (4.4), 
(4.5), both Si,n and S'2,n, can be expressed as a sum of 2n terms each having the expo- 
nential structure ± exp{a + bw} lL4-<ui<£/t for appropriate constants a, b varying among 
the 2n terms. Thus, the quantity in (5.4) can be expressed as: 

fci.„ 

1 + ^{-If' exp{a., + lL4-<^<(7t 
1=1 

for ki^n — 4{(2n + 1)^ + (2n + 2)}, and constants a^, &i, Ci with Ci G {0, 1}. Working 
similarly for all four summands forming up 5'^„+i in (5.2) we get that the function 
fin+iiw) can in fact be written as the weighted sum: 

/2n+i(w) = '^{-'^y exp{ai + biw}lL,<w<Ui x 7r(w) (5.5) 

i=i 

for fc„ = 2(fci^„ -|-fc2.n) with fc2,,i = 4{(2n-|-2)^-|-(2n-|-l)}, and some explicit constants a^, 
bi,Ci S {0, 1}, Li, Ui. Experimentation has showed that /i is already a very good envelope 
function for the rejection sampler, in which case fc„ = fco = 64; this is not accidental, 
and relates with the rapid exponential convergence of the alternating sequence in (4.5) 
to its limit. The cdf, say Fi(w), corresponding to the unormalised density function fi{w) 
can be analytically identified since integrals for each of the summands in (5.5) can be 
expressed as differences of the cdf of the standard Gaussian distribution. Samples from 
/i(u;) can then be generated using the inverse cdf method, i.e. by returning F^^{R) 
for R ^ Un[0, 1]. cannot be found analytically, but numerical methods can return 
F^^{R), up to maximum allowed computer accuracy, exponentially fast. We have used 
MATHEMATICA to automatically calculate all integrals giving the cdf, and then incorporated 
the calculation into a C++ code. 

Summarising, our rejection sampler will be as described below, where for simplicity 
we write p{'w) = p{L^, L'^, U^, U^; q, r, x, w, y): 
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Repeat until the first accepted draw: 

Propose w ^ fi and accept with probability f{w)/fi{w) = p{w) / Sf (w) . 

Note here that the acceptance probability involves p^w) which is made up of eight infinite 
series, see (4.9). We avoid approximations by using the alternating construction (5.3) 
and employ the methods of Section 4.2 to obtain the value of the decision indicator 
^R<piw)/sfiw) for some i? ~ Un [0, 1]. 

As shown in Step 1 of Table 1, once Xf is obtained, we adjust the allowed range for 
the extrema of the bridge on [s, t] by simply setting ~ Xf , = A Xt* . 

5.2. Bisection of Xgf: Updating the Layers Given Xf. 

At the second step of the bisection procedure, see Table 1, we obtain separate information 
for the extrema of the two newly formed bridges given the middle point Xt* : the one 
bridge being from Xg to Xf, the other from Xf to Xt. In particular, the algorithm 
will decide over the range of the four newly formed layers, Cg.f, Us.f, Cf^t, l^t*,t in the 
following manner: for the case of Cs^f for instance a decision will be made over whether 
nis^f lies in [L^ , Li^] (which is the allowed range for the minimum of the original bridge 
on [s, t\) or in [L^, Xs hXf]. The apparent analogues apply in the case of the three other 
layers. 

One might initially think that there are in total 2'^ different scenaria for the four 
layers. But one has to remember that the update has to respect the information in X^^t, 
so that at least one of the two minima (resp. maxima) on [s,t*] and must lie in 

[L^^L^] (resp. \U^,U'^]). In particular, there are in fact nine different possible scenaria, 
which are the ones shown in Table 3 (labelled as events {E = i}, for 1 < i < 9): a 
value of 1 in Table 3 means that the corresponding minimum or maximum will still 
be found within the allowed range for the original bridge on [s,t], whereas a value of 
means that the second option occurs and the extremum will be shifted inwards. For 
instance, a value of for the indicator variable concerning TOs,t*, M^.f, ruf^t or Mf^t 
implies that TO,,t. e [L^,Xs^Xf], Ms,f € [X^V Xf,U^], vfif^t S [L^,Xf A Xt] or 
Mf^t e [Xf y Xt,U^] respectively. 

The probability for each of the events in Table 3 can be derived via functions /3(-; •) 
and p{-; •) defined in (4.8) and (4.9) respectively. Recall that we are conditioning upon 
Xs^t and Xf , so we work as follows: 

V[E = i\I,^uXf]^V[E = i\m(l [L\ L\ M e U^, X,, Xf , X^ ] 
^ P[E = i,m€ [L\ L\M& [U^, U^] \ X„ Xf , X^ ] 

P [m e [L^Lt], M G I A„ Xf, Xt] 

^ V[E = i\ Xs,Xf,Xt] 
p(L^,L'^,U^,U^;q,r,x,w,y) 

Now, conditionally on {Xs, Xf , Xt} the law of the path factorises into two independent 
Brownian bridges. Thus, recalling also the definition of /?(•; ■) in (4.8), the probability 
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~P[E = i \ Xs, Xf , ] in the numerator above can be written as a product of two /?(•; •) 
functions. The analytical calculation of the numerator, or equivalently of the product 
¥[£ = i\ Is^t,Xf ] X p{L^, i^, U^, f/^; q, r, x, w, y), is given in Table 4 where, to simplify 
the presentation, we have set: 

Wx = X /\ w, = xV w, Wy ^ w A y, ^ wW y . 

The method to simulate the discrete random variable E could follow the alternating 
series approach of Section 4.2. Analytically, consider the cumulative probability values 
Pi = F [E < i \ Xs^t,Xt*]. A simple inverse cdf sampling method requires finding the index 
inf{i > 1 : i? < pi} for a i? ~ Un [0, 1]. Note now that the pi's can be bounded above and 
below by monotone converging sequences as in (4.10), thus each comparison {R < pi} 
can be carried out via the alternating series approach of Section 4.2 without any need 
for approximations. 



Event 


Left bridge 




Right bridge 




E = i 




t,eiui-,ut] 




e[(7J-,(7t] 


i = 1 


1 


1 


1 


1 


i = 2 


1 


1 





1 


i = 3 


1 


1 


1 





i = 4 


1 


1 








i = 5 





1 


1 


1 


i = 6 





1 


1 





i = 7 


1 





1 


1 


i = 8 


1 








1 


i = 9 








1 


1 



Table 3. The nine possible scenaria for the extrema of the two Brownian bridges (from Xs to X* and 

from Xt* to Xt). 



i 


PlE = i\Is,t,X 


t* 




X, 


w,y) 


1 


I3{L-^,L^,U^,U^ 


1 


x,y) 


X /3{L^,LT,U^,UT 


r, 


w,y) 


2 




q, 


X, w) 


X l3(L^,Wy,U-^,U'f 


r 


w,y) 


3 




q, 


X, w) 


X l3{L^,L'^,wy,U^ 


r 


w,y) 


4 


I3{L^,L^,U^,U^ 


q, 


X, w) 


X l3{L'^,Wy,wy,U^ 


r 


w,y) 


5 


l3iL-^,w^,U\Uf 


q 


X, w) 


X p(L^,L'f,U-^,U^ 


r 


■w,y) 


6 


I3{L^,L'^,U^,U^ 


q, 


X, w) 


X l3(L-^,L'f,U-^,U'f 


r 


w,y) 


7 




q 


X, w) 


X l3{L^,L^,wy,U-^ 


r 


■w,y) 


8 




q, 


X, w) 


X l3(L^,Wy,U^,U'f 


r 


w,y) 


9 




q 


X, w) 


X I3{L^,L'^,U-^,U^ 


r 


w,y) 



Table 4. The conditional probabilities for each of the events in Table 3. 



5.3. Remaining Sampling Procedures 

A sampling algorithm is required for the refinement of the uncertainty over the extrema 
of a Brownian bridge. As described in Section 2.1, given the current intersection layer 
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information Ts,t and in particular the fact that M^.t € [L^, L^], the algorithm will need to 
decide whether the maximum M^^t lies in [U^, U*] or in [U*,U^], for U* = {U^ + U^)/2. 
Recalling the definition of /3(-; •) from (4.8), it is easy to check that the ratio: 



provides precisely the probability of the event {Mg.t G [U* , J7^] |Is.t}. Thus, wc can again 
use the alternating sequence construction of Section 4.2 to simulate, without approxima- 
tion, the binary variable ^ m^, t&[u* .uf]- The same approach can be followed for refining 
the allowed range for the minimum TOs,*- 

We should also give some details over the initialization of the layers Z^o.i and £o,i at 
the first step of the e-strong algorithm in Table 2 given Xo and Xi. (Note that sometimes, 
as in the example applications that we consider in the following section, this initialization 
steps might not even be necessary, as the problem at hand provides a natural definition of 
Uq^i and Cq^i.) One way to proceed is by specifying increasing sequences {ai}i>Q, {6i}i>o, 
with ao = 6o = 0, growing to oo and a bivariate index / such that: 

{/ («, j)} = {x -Oi <m<x - a^-i, y + bj-i < M <y + bj} , 

where x = x/\y,y = x\/y. Wc can easily identify the probability distribution of / under 
the Brownian bridge dynamics since: 



Thus, we can work as in the case of the simulation of the discrete variable E in Section 5.2: 
assuming / = 1,2,... is some chosen ordering of the states of /, an inverse cdf method 
would required finding inf{i >l:i?<P[/<z]} for R Un[0, 1], and approximations 
at the comparison between R and P [ / < i] can be avoided via the alternating series 
approach. In practice, one could select some big enough values for the first couple of 
elements of the sequences {oi} and {hi] so that almost all probability mass is concentrated 
on {/ = (i, j)} for i, j < 2, and not a lot of computational resources are spent on this 
step. 

6. Application: 

Unbiased Estimation of Path Expectations 

The information provided by the e-strong algorithm can be exploited to deliver unbiased 
estimators for path expectations arising in applications, avoiding discretization errors 
characterising standard approaches. We emphasize that we mean to sketch here only a 
potential direction for application of the algorithm. Analytically, consider a non-negative 
path functional F : C([0, 1],M) ^ and the expectation: E [F{X)], X being a Brow- 
nian motion on [0, 1]. One can easily check, by integrating out that: 



P{L\L\U\U^-l,x,y) 
P{L^,L^,U^,U^-l,x,y) 



x -ai,x - aj_i, y + fej-i, y + bj; l,x, y) . 



^F(X)>E • e 



E ^ Exp(l) , 



(6.1) 
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with E being independent of X, is an unbiased estimator of E[F(X)]. The e-strong 
algorithm could be utilised here to unbiasedly obtain the value of the binary variable 
I F(x)>E in finite computations. We can easily find the second moment of the unbiased 
estimator in (6.1): 

E[If(x)>b -e^^l-Eie^f^'l-l . (6.2) 

We describe for a moment in more detail the identification of ^f(x)>u ^'^^ the e-strong 
algorithm. Utilising the lower and upper convergent processes X'^{n), X^(n) in (3.1) 
one could in many cases analytically identify quantities F^, (realisable with finite 
computations) such that: 

Ft < Ft+, < F{X) < Fl^, < Ft ; 
Ft - Fi ^ . 

Given enough iterations, there will be agreement; for the a.s. finite random instance: 

K = inf{n > : = I^t>£; } (6.3) 

we will have 

Thus, combining (6.1) with (6.4), we have developed an unbiased estimator of a path 
expectation, involving finite computations. Certainly, the numerical efficiency of such an 
estimation will rely heavily on the stochastic properties of k and the cost of generating 
Fyj;, Ft, and of course the variance of the estimator. 

The particular derivation of the above unbiased estimator of the path expectation is 
by no means restrictive; one can generate unbiased estimators using distributions other 
than the exponential. Consider the following scenario. We can generate some preliminary 
bounds Ft^, Ft^^ up to some fixed or random (depending on X) instance uq. Now, one 
can easily check (by considering the conditional expectation w.r.t. R\X) that: 

lFix)>R Fl + iFixxRFi ; R ^ \Jn[Fi , Ftj , (6.5) 

is also an unbiased estimator of E [ F{X) ] . We have empirically found the estimator (6.5) 
to be much more robust than (6.1) in the numerical applications wc present in the sequel. 
This is not accidental: for instance, considering a random no such that Ft^ — F^^ < C, 
for a constant C > 0, we get that the second moment of the estimator (6.5) will be: 

^[lFiX)>R {Flf + lFiX)<R{Fif ] = 

= E [ Fix) (Fl + Ft J ] - E [ Ft^Ft^ ] < E [ F^ (X) ] + C E [ F(X) ] 

which has now a quadratic structure - compare this with (6.2). In general, increasing 
no adds to the computational cost per sample, but decreases the variance. We have 
empirically found that moderate values of no deliver significantly better estimates than 
(6.1) and will be using such an approach for our numerical examples in the sequel. 
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6.1. Numerical Illustrations 

We will apply the e-strong algorithm to unbiasedly estimate some option prices arising 
in finance. In particular option prices arc expressed as expectations: 

E[F(5)] 

of a functional F{-) of the path process S = {St} modelling the underlying asset. We will 
consider some double-barrier options corresponding the expectations of the functionals: 

Fa{S) = e"'~'^(supS't - Ks)^ '^Ls<infSt<supSt<Us (6-6) 

Fb{S) = (y / Stdt-Ks) I Ls <inf St <supSt<c/s ; (6-7) 
Jo 

Fc(S') = e"''^(supS't - A's)+ lLs<i„fSte-'-*<supSte-'-*<C/s ' (6-8) 

(where for x € R, x V 0) for underlying asset S = {St;t £ [0,r]} modelled via a 

geometric Brownian motion (we consider a Black and Scholes framework) determined as: 

log St ^log So + {r-2^)t + aWt (6.9) 

for constants r (interest rate), cr > (volatility) and a Brownian motion {Wt}. Also, T 
above is the maturity time, Ks the strike price and Ls, Us the lower and upper barriers 
respectively; suprema and infima are considered over the time period [0,r]. Note that 
£[^6(5*)] corresponds to the price of the Asian option, see e.g. [21]. The process St is 
an 1-1 transformation of a Brownian motion with drift. In particular, we can rewrite the 
functionals (6.6)-(6.8) as follows: 

i^,(X) = e-''^(e'^™P^' -1^5)+ lL<infx,<supX,<c/ ; (6.10) 

FbiX) = e-^^{ i / e"^' dt - Ks)^ lL<inix,<snpX,<u ; (6.11) 
Jo 

F,{X) = e-'^^(e'^^^P(^*+^') -/is)+ lL<i„f x,<«upX,<c/ , (6.12) 

for L = log(Ls)/(T, U = log([/s)/cr, and: 

Case Fa, Ft : Xt = log(S'o)/a + (^ - §)t + Wt ; 
Case F,: Xt = log(S'o)/a -^t + Wf 

Conditionally on its ending point, the dynamics of the drifted Brownian motion do not 
depend on the value of the drift and coincide with those of a simple Brownian bridge; 
this is a simple by-product of the Girsanov theorem, see e.g. [15]. Thus, the e-strong 
algorithm can now deliver convergent, lower and upper dominating processes for X. 

The choices of functionals in (6.6)-(6.8) is not accidental: some generic characteristics 
of the structure of each functional (relevant also for other applications) will effect the 
set-up of the e-strong algorithm and its efficiency; we will say more on this in the sequel. 
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For all three examples, our general methodology is as follows: we begin by sampling Xt 
and, then, the indicator variable lL<inf Xt<supXt<c/; if the latter is the sample for our 
unbiased estimator is simply 0, otherwise we proceed with applying the methods of the e- 
strong algorithm by initializing the first intersection layer as Iq.t = {-'^o, XtMo.t, -^ct} 
for intervals Wo,t = [-'^o V Xt,U] and £o,t = [L^Xo A Xt\- In some cases we might 
not need all of the machinery of the e-strong algorithm to construct the sequences F^^ 
enveloping F{X), with direct implications on the efficiency of the algorithm, as we 
explain analytically below. 

Fa-example: Only Refinement 

Here, we need information only on the marginal variable supXj (and not the whole of 
the continuous path on [0, T]) to develop an alternating series for Fa{X). Thus, it suffices 
to apply a reduced version of the complete e-strong algorithm in Table 2 where we only 
repeatedly refine the initial intersection layer Xq.t (in particular, we only refine the layer 
for the maximum) as described in Section 2.1 (and never bisect it) to construct F^, F^. 
In particular, having defined: 

q^{x)=e--^ie-- -Ks)+ . 

knowing that after n refinements the allowed range for supXt is [U^,UX\ (with initial 
position [U^, U^] = [Xo V Xt, U]) wc set F^ = (biU^), F^ = (/.(f/t). 

Fh-example: Refinement And Bisection 

The complete machinery of the e-strong algorithm in Table 2 is required here as we need 
to bound a path integral. Recall that the n-th step of the algorithm provides the piecewise 
constant paths xf{n), xj{n) enveloping X defined in (3.1). We now set F^ = Fb{X-^{n)), 
Fl^FMHn)). 

Fc-example: Selective Refinement And Bisection 

Wc will now only need to bisect a selection of intersection layers as wc will be allowed to 
delete intersection layers that cannot definitely contain sup(^ t+Xt) during the execution 
of the e-strong algorithm. In particular, assuming the current collection (after n — 1 steps) 
of stored intersection layers Tsi.u, with Si <ti < s,;+i, containing information about the 
path X, and determining the allowed range for sup(^ t + Xt) : 

[Ui-^.UI_^] = [^MUtu + ^MUl,u + 7 (6-13) 

i i 

wc proceed to the n-th step where: (i) wc bisect and refine all stored intersection layers 
Isi.ti, (ii) calculate the running bounds [t/^, U^] by taking the suprema as in (6.13) but 
now over all newly obtained intersection layers, (iii) delete the obtained intersection layers 
Is,t for which 11}^ < (as they cannot offer extra information on the whereabouts of 
sup(^ t+Xt) given that we already know that the latter is within [ U^, Uj^]) and store only 
the remaining ones for the next iteration. At each step we set F^ = 4){U}[), F^ = 4>{UX) 
with (j) as defined above. 
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Numerics 

We have run the e-strong algorithm for the above scenaria. To give an idea of its comput- 
ing cost, we compare its execution times with those of the standard (Euler) approximation 
method that replaces the continuous-time path {St t € [0,T]} with its discrctiscd ap- 
proximation {Sti}\^Q, for ti = + 5, with step-size S = T/l; then, continuous-time 
maxima and integrals appearing in the functionals (6.6)- (6.8) are replaced with their 
obvious approximations based on the discrete-time vector {S'ti}- We run our simulations 
under the parameter selections: 



0.05 



0.2,Sa = l,K = l,T = l,U = 1.25 , L ^ 0.75 



Tables 5-7 show results from the simulation study. For each different algorithm we show 
its execution time (all algorithms were coded in C++) and a 95% confidence interval for 
the mean of the realised estimators to give an idea about the variance of the estimates 
and their bias (for the case of the Euler approximation, as the e-strong algorithm is 
unbiased). The results in Tables 5, 7 are obtained via 100,000 independent realizations 
of the estimators, whereas those in Table 6 via 10,000 independent realizations. 



Euler Approximation 


<5 


time (sees) 


95% Conf. Int. 


1/10 


0.4 


[638, 647] X 10— ' 


1/20 


0.8 


[657, 667] X 10""' 


1/40 


1.5 


[669, 679] X 10""' 


1/80 


2.9 


[674, 683] X 10-4 


1/160 


5.8 


[680, 689] X 10-'' 



e-strong 


no 


time (sees) 


95% Conf. Int. 


2 


1.1 


[683,693] X 10-" 



Table 5. Simulation results from the application of the Euler approximation and the e-strong 
algorithm for the estimation of the option price in E [Fa{S) ] in (6.6). The results in the table 
correspond to a sample of 100, 000 estimates. 5 is the discretisation increment of the Euler method, and 
no is the number of preliminary steps for the e-strong algorithm before the simulation of the uniform 

random variable (see (6.5)). 



Euler Approximation 


S 


time (sees) 


95% Conf. Int. 




0.1 


[157,169] X 10-" 


10-2 


0.4 


[120,130] X 10-" 


10-3 


3.5 


[106,116] X 10-" 


10-" 


34.3 


[107,116] X 10-" 


10-5 


344.8 


[102,112] X 10-" 



e-strong 


no 


time (sees) 


95% Conf. Int. 


2 


115.9 


[81,128] X 10-" 



Table 6. Similar results as for Table 5, but now for the ease of the Asian option E [ -F(,(5) ] in (6.7) 
with the difference that here the results correspond to a sample of 10, 000 estimates. 



Looking at the three tables, we can make some comments; we focus more here on 
giving a simple picture to the reader than being mathematically precise. The cost per 
sample of the e-strong algorithm corresponds to that of the Euler approximation with 
5 1/40, 5 « 4-1 • 10-" and 5 « 10-3 for the cases of E [ (£■) ] , E [ i^,, (S*) ] and E [ (5) ] 
respectively. Taking also the standard deviation under consideration (but not the bias) 
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Eulcr Approximation 


<5 


time (sees) 


95% Conf. Int. 


1/10 


0.5 


[797, 807] X 10-" 


1/20 


0.9 


[822,832] X 10-"' 


1/40 


1.7 


[833, 844] X 10"'' 


1/80 


3.3 


[835, 846] X 10"'' 


1/160 


6.6 


[846, 858] X 10"'' 



e-strong 


no 


time (sees) 


95% Conf. Int. 


2 


178.0 


[842,854] X 10-" 



Table 7. Similar results as for Table 5, but now for the case of E [i^c{S) ] in (6.7). 



from the column with the confidence intervals, for the case of E \ Fb{S) ] we would need 
about 25 times more samples than then Euler approximation to attain the same range 
for the confidence interval; thus, ignoring the bias for the Euler approach, one could say 
that the overall cost of the e-strong algorithm for the case of E [-F'fc(S') ] corresponds to 
that of the Euler method with step-size J' w (4 • 2b)~^lQ~^ = 10"^. 

However, a general remark here is that the e-strong algorithm returns unbiased estima- 
tors of the relevant path expectations, and for the applications we have considered above 
it can provide accurate, unbiased estimates in reasonable amounts of time. Even when 
ignoring the bias of the Euler approach, for the cases of E [ Fa{S) ] and E [ Fc{S) ] the cost 
of the e-strong algorithm already seems to be on a par with that of the Eulcr method for 
relatively non-conservative choices of discretisation step 5. (We should also stress that 
there is definitely great space for improving the efficiency of the used computing code for 
the e-strong algorithm.) 



6.2. Remark on Number of Bisections for e-Strong Algorithm. 

We make a comment here on the number of required iterations before the value of the 
binary variable I f{x)>r in (6.5) is decided. Proposition 3.2 will be of relevance in this 
context. Recall that k in (6.3) denotes the number of steps to decide about ^f{x)>r- 
The cost of K iterations of the e-strong algorithm (when its full machinery is required) is 
proportional to K. = 2'^ . In the context of (6.5), we find: 

P[/C>2"]^E[P[K>n|X]]^E[ ^""^"1 ] . 

Proposition 3.2 states that \X'^{n) — X^(n)\L^ = 0(2""/^). The same rate of conver- 
gence will many times also be true for E [ — ] : this will be the case for instance 
when F{X) = /(/^ g{Xs)ds) under general assumptions on f,g (e.g. if \f{y) — f{x)\ < 
M{x,y)\y — x\, for a polynomial M, and the same for g; a proof is not essential here). 
For such a rate (and since the user-specified F^^ ~ F^j^^ should be easily controlled) we 
will get: 

P [/C > 2"] = 0(2-"/2) 
giving the infinite expectation E [/C] = oo. 
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In a given application though, one could fix a big enough maximum number Umax, 
stop the bisections if that number has been reached and report, say, {F^^ + F^^)/2 as 
the realization of the estimator if that happens, without practical effect on the results. 
To explain this, note that we know, from (6.5), that the actual unbiased value is either 
or F^^ , so we know precisely that the absolute bias from the single realization when 
nmax was reached cannot be greater than [F^^ — F^^)/2. In total, when averaging over 
a number of realizations we can have a precise arithmetic bound on the absolute value 
of the bias of the reported average; if Umax is 'big enough' so that we reach rimax only 
in a small proportion of realizations the (analytically known) bias could be of such a 
magnitude that the reported results will be precisely the same as when implementing 
the regular algorithm without rimax for a reasonably selected degree of accuracy. For 
example, in the case of the estimation of ¥j[Fc{S)] in Table 6, we have in fact used 
nmax = 10 and found that the introduced bias was smaller than 3 x 10^^ so avoiding it 
would not make any difference or whatsoever at the results reported right now in Table 6. 

Note that such an issue did not arise in the cases of E [ Fa{S) ] and E [ Fc{S) ] when a 
reduced version of the e-strong algorithm was applied. 

7. Further Directions for Applications 

We sketch here some other potential applications of the e-strong algorithm. 

In the case of barrier options, sometimes one needs to evaluate expectations involving 
a Brownian hitting time (see e.g. [19]). Given a non-constant boundary H : [0, oo) — >■ R, 
such that So < Hq, consider: 

T„ = M{t >0:St>Ht}, 

with S = {St} being the geometric BM in (6.9). The price of a related derivative will be 
E [ F{S) ] where now: 

F(5)=V(5t)-I.^<t 

for some pay-off function ip{-). This estimator requires the evaluation of It^<t for a 
realised path. Such an evaluation is possible under our simulation methods, since for a 
given bridge, say from Ss to St with s < t < T, we can decide if its maximum is within 
[H^^,Hlf] or not (thus, deciding also whether there is a chance that the bridge hits H 
on [s,t] or not), with 

= mi{Hu; u G [s, t]}, = sup{F„; u e [s, t]} , 

using the refinement procedure described in Section 2.1 (more particularly, a slightly 
modified version of it, where instead of halving the allowed variation for the maximum, 
it decides if it lies in a given interval or not). Computational effort will then only be spent 
on the bridges for which the maximum is indeed in [H^ f., iteratively bisecting them 

until a definite decision is reached about whether H has been hit. 
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Individual simulation techniques employed in the development of the e-strong algo- 
rithm are also of independent interest. For instance, we have exploited during the con- 
struction of the e-strong algorithm a monotonic property at the core of the Brownian 
structure; we can further use this characteristic to develop original simulation techniques 
for Brownian distributions. One application for instance could involve dynamics of Brow- 
nian motion restricted to stay in a bounded domain. A Brownian motion with constant 
drift, restricted to remain in (— 7r/2, 7r/2), is known (see [16]) to be described via the 
stochastic differential equation; 

dXt = - lmi{Xt)dt + dWt . (7.1) 

Unbiased sampling methods for Xt are not (to the best of our knowledge) available; one 
has to resort to Euler, or other, approximations. We can, however, now construct an 
exact sampling algorithm based on the methods so far described. Girsanov's theorem 
provides the following expression for the transition density of the Markov process (7.1): 

p{y;x,t) ■. = V[Xtedy\ Xq = x] / dy 

= )^-f{-7^/2,^^/2;t,x,y)po{y■,x,t), y G (-7r/2, ^2) , (7.2) 

cos(a;) 

for the unconditional Brownian transition density: 

Po{y;x,t)^{2nt)-'/'e-^y-^^'/(^'^ . 

Density (7.2) has a structure reminiscent of that of the density of the middle point in 
Proposition 5.1: ideas employed there, are also relevant now. Analytically, for x not close 
to the boundaries, one can simply carry out a rejection sampler with proposals from po- 
Then, the acceptance/rejection decision will be based on comparing a real number with 
7(— 7r/2, 7r/2; <, x, J/) following the pattern described in Section 4.2. As x approaches the 
boundaries, this algorithm becomes inefficient. But, similarly to the method for the sim- 
ulation of the density in Proposition 5.1, partial sums from the infinite series-expression 
for 7(— 7r/2, 7r/2; t, x, y) can be incorporated in the proposal to produce an efficient algo- 
rithm (Section 5.1 describes the algorithm for the more complex density appearing there; 
here, we omit the details). 



8. Conclusions 

We have presented a contribution to sampling methods for Brownian dynamics: a new 
iterative algorithm that envelopes the Brownian path, thus offering explicit information 
for all its aspects (minimum, maximum, hitting times). Individual steps of the algorithm 
could be of independent interest, yielding new sampling methods for distributions derived 
from Brownian motion dynamics. 

We should remark here on the generality of our scope. The e-strong algorithm (or 
some of its individual steps) can provide, more or less unchanged, unbiased Monte-Carlo 
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estimators in separate estimation problems, for which quite an extensive amount of case- 
specific methods have been investigated in the hterature; for instance, one can refer to 
the long literature for the applications we briefly described in Sections 6 and 7. 

We have presented some applications and sketched some others towards illustrating 
the potential of our methods. Note that the infinite expectation issue remarked in Sec- 
tion 6.2 is a direct consequence of the Brownian dynamics: the maximum of the Brownian 
path scales as At^/^ on a small time interval [0, At] (see e.g. [18]). Thus, any enfolding 
processes will necessarily converge not faster than 0(At^/^) (which is the order attained 
by the e-strong algorithm). This relatively slow convergence of the enfolding processes 
also explains the increased cost for when estimating E [Fh{S) ] in Section 6.1. We envis- 
age that it might be possible to combine the iterative process of the algorithm with a 
coupling step once the bounding processes are relatively close to each other to overcome 
this long anticipation (in the spirit of [20], [4]). We hope to formalise this idea in future 
research. 
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